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1/^ ' 

A practical high-accuracy relativistic method of atomic structure calculations for univalent atoms 
| is presented. The method is rooted in the coupled-cluster formalism and includes non-perturbative 

^v^j , treatment of single and double excitations from the core and single, double and triple excitations 

involving valence electron. Triple excitations of core electrons are included in the fourth-order of 
many-body perturbation theory. In addition, contributions from the disconnected excitations are 
incorporated. Evaluation of matrix elements includes all-order dressing of lines and vertices of the 
diagrams. The resulting formalism for matrix elements is complete through the fourth order and 
sums certain chains of diagrams to all orders. With the developed method we compute removal 
energies, magnetic-dipole hyperfine-structure constants A and electric-dipole amplitudes. We find 
i | . that the removal energies are reproduced within 0.01-0.03% and the hyperfine constants of the 3si/2 

and 3pi/2 states with a better than 0.1% accuracy. The computed dipole amplitudes for the principal 
Oh 3si/2 — 3pi/2 ; 3/2 transitions are in an agreement with 0.05%-accurate experimental data. 

2 PACS numbers: 31.15.Dv, 31.30.Jv, 32.10.Fn, 32.10.Hq, 32.70.Cs 

o : 

I. INTRODUCTION 

C/2 ■ 
U ■ 

This work is aimed at designing a practical ab initio atomic-structure method capable of reaching accuracy at the 
^-y level of 0.1% for properties of heavy univalent many-electron atomic systems. The improved accuracy is required, for 
i-C ' example, for a refined interpretation of atomic parity violation (APV) with atomic Cs and planned experiment 

with Ba + . At present namely the accuracy of solving the basic correlation problem is the limiting factor in the 
APV probe of "new physics" beyond the standard model of elementary particles. In addition, it is anticipated that 
t-H , the improved accuracy would unmask so far untested contributions from quantum electrodynamics (QED) in heavy 
^ • neutral many-electron systems 0. 

Here we report developing a many-body approach based on the coupled-cluster (CC) formalism In the CC 

formalism the many-body contributions to wave function are lumped into a hierarchy of multiple (single, double, ...) 
particle-hole excitations from the lowest-order state. Due to a computational complexity, previous relativistic CC-type 
calculations [1,0,0,0,0,01 for univalent atoms were limited to single- and double excitations. Triple excitations 
' were treated only in an approximate semi-perturbative fashion 0, 0, II^ . 0, llil Hlj . Compared to these previous 
calculations, here we fully include valence triple excitations in the CC formulation; we will designate our approximation 
""^5 ' as CCSDvT method. Further, compared to calculations by Notre Dame group, here we also incorporate a subset of 
so-called disconnected excitations (non-linear CC terms). For sodium atom, such non- linear CC terms were previously 
' 55 ' included in Ref. |Toj and in non-relativistic calculations |16| . Finally, in calculations of matrix elements we include 
CC-dressing of lines and vertices [l7j and we also directly compute complementary fourth-order diagrams (mainly 
. due to core triple excitations). The resulting formalism for matrix elements is complete through the fourth-order of 
many-body perturbation theory (MBPT) and also subsumes certain chains of diagrams to all orders. 
;> As a first application of our method, we carry out numerical calculations for atom of sodium. Sodium (11 elec- 

trons) has an electronic structure similar to cesium (55 electrons), but it is not as demanding computationally. By 
computing properties of Na atom we observe that a simultaneous treatment of triple and disconnected quadruple 
■ excitations is important for improving theoretical accuracy, as the two effects tend to partially cancel each other. We 
compute removal energies, magnetic-dipole hyperfine-structure (HFS) constants A and electric-dipole amplitudes for 
the principal 3si/2 — "&Pj transitions. We find that the removal energies are reproduced within 0.01-0.03% and the 
HFS constants of the 3s and 3pi/ 2 states with a better than 0.1% accuracy. The computed dipole amplitudes are 
in a perfect agreement with the 0.05%-accurate experimental data. However, our result for the HFS constant of the 
3^3/2 state disagrees with the most accurate experimental values [l8L fli^ by 1%, while agreeing with less accurate 
measurements [201121). 

The paper is organized as follows. First we discuss generalities of the coupled-cluster formalism and many-body 
perturbation theory in Sec. [H] Explicit CCSDvT equations and analytical expressions for energies, matrix elements, 
and normalization corrections are presented in Section lTTTl In Sec. II VI we tabulate and analyze the results of numerical 
calculations of properties of sodium atom. Finally, we draw conclusions in Sec. Unless specified otherwise, atomic 
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units |e| — h = m e — Atteq = 1 are used throughout. 

II. GENERALITIES 

In this Section we recapitulate relevant formulas and ideas of atomic many-body perturbation theory (MBPT) and 
the coupled-cluster formalism for systems with one valence electron outside the closed-shell core. 

A. Atomic Hamiltonian and conventions 

The Hamiltonian of an atomic system may be represented as 



# = ( E W r *) + E ^HF(ri) ) + I \ E J 



— -E^hfW , (1) 

ri l i J 

where h nuc is the Dirac Hamiltonian including kinetic energy of electron and its interaction with the nucleus, £7dhf 
is the Dirac-Hartree-Fock (DHF) potential, and the last term represents the residual Coulomb interaction between 
electrons. To reduce the number of MBPT diagrams, we employ frozen-core (or V^ -1 ) DHF potential j^. The 
single-particle orbitals ipi and energies £j are found from the set of DHF equations, 

(ftnuc + ^DHF) <fi = EUfi . (2) 

The Hamiltonian in the second quantization reads (omitting common energy offset) 

H = H + G = Y^ e<JV[otoi] + ^ E 9im N \4 a \ a l a k] , (3) 

i ijkl 

where operators a, and a\ are annihilation and creation operators, and N[- ■ ■ ] stands for a normal product of operators 
with respect to the core quasi- vacuum state |0 C ). Labels k and I range over all possible single-particle orbitals. In 
the following we will employ a labeling convention where letters a, b, c are reserved for core orbitals, indices m, n, r, s 
label virtual states, and letters v and w designate valence orbitals. In this convention valence orbitals are classified 
as the virtual orbitals. In Eq. the quantities gijki are two-body Coulomb matrix elements 

Qijkl = Jd 3 rJ d3 rV t (r) ^y)_L ^(r') . (4) 

Notice the absence of the one-body contribution of G in the second-quantized Hamiltonian, Eq. ; this simplifying 
feature is due to the employed V 1 ^^ 1 approximation and leads to a greatly reduced number of terms in the CC 
equations. 

In MBPT the first part of the Hamiltonian © is treated as the lowest-order Hamiltonian Hq and the residual 
Coulomb interaction G as a perturbation. In the lowest order the atomic wave function with the valence electron in 

an orbital v reads = |0 C ) . Further the wave operator £1 is introduced; it promotes this lowest-order state to 

the exact many-body wave function 

|^>=fi|*(°>>. (5) 

In the conventional order-by-order MBPT, a perturbative expansion for operator f2 is built in powers of residual 
interaction G resulting in a hierarchy of approximations for correlated energies and wave-functions. 

B. Coupled-cluster method 

One of the mainstays of practical application of MBPT is an assumption of convergence of series in powers of the 
perturbing interaction. Sometimes the convergence is poor and then one sums certain classes of diagrams to "all 
orders" using iterative techniques. The coupled-cluster formalism is one of the most popular all-order methods. The 
key point of the CC method is the introduction of an exponential ansatz for the wave operator pflj 

Q = N[exp{K)] = 1 + K + ^N[K 2 } + . . . , (6) 
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where the cluster operator K is expressed in terms of connected diagrams of the wave operator f2. The operator K is 
naturally broken into cluster operators (K) combining n simultaneous excitations of core and valence electrons from 

the reference state |*£ 0) ) to all orders of MBPT 



total number of electrons 

E 



K= V {K)=S + D + T + ..-, (7) 



i.e., K is separated into singles (S = (K)^, doubles (D = (K) 2 ), triples (T = (K) 3 ), etc. For the univalent systems 
we further separate the cluster operators into two, core and valence, classes 

(K) n = (K c ) n + {K v ) n . (8) 

Clusters (K c ) n involve excitation from the core orbitals only, while (K v ) n describe simultaneous excitations of the 
core and valence electrons. Then S = S c + S v , D = D c + D v , etc. 

A set of coupled equations for the cluster operators {K) n may be found from the Bloch equation specialized 
for univalent systems [24| 

(s v -H )(K c ) n = {QGn} conncctcdn , 
(e v + 5E V - H ) (K v ) n = {QG ^} conncctcdi „ , (9) 

where the valence correlation energy 

8E V = (¥ v V\Gn\¥V), (10) 

andQ = l-|*i 0) )(^ 0) | is a projection operator. Notice that only connected diagrams are retained on the r.h.s of the 
equation, r.h.s. diagrams being of the the same topological structure as clusters (K) . The resulting CC equations 
for the core clusters do not depend on the valence state. 

Although the CC approach is strictly exact, in practical applications the full cluster operator K is truncated at a 
certain level of excitations, e.g., at single and double excitations (CCSD method). In particular, for univalent atoms 
the CCSD parametrization may be represented as 

K SD = S c + D c + S v + D v = 

ma mnab m^v mna 

The cluster amplitudes p... are to be determined from the Eq.©. 

A linearized version of the CCSD method discards non-linear terms in the expansion of exponent in Eq. © of the 
coupled-cluster parametrization, i.e., Q SB = 1 + K su . This leads to discarding disconnected excitations from the 
exact many-body wave function. We will refer to this approximation simply as singles-doubles (SD) method. For 
alkali-metal atoms the SD method was employed previously by the Notre Dame group H, |j| 0, . The resulting SD 
equations are written out in Ref . ■ A typical ab initio accuracy attained for properties of heavy alkali-metal atoms 
is at the level of 1%. 

Successive iterations of the CC equations I© recover the traditional order-by-order MBPT. As discussed in Ref.[j|, 
the core and valence doubles appear already in the first order in the residual interaction G: 

^ 9mnab 

Pmnab ~ j , \- ) 

&a \ &b &m &71 
^ 9mnva 

Pmnva ~ j ■ \ J 

£v i Sa £m £n 

Valence and core singles appear at the second iteration of the CC equations and are effectively of the second order in 
G. We will employ this "effective order" classification to develop our approximation to the CC equations. 



C. Triple excitations. Motivating discussion 



Certainly the truncation of the CC expansion leads to a neglect of many-body diagrams containing excitations 
beyond singles and doubles. For example, both the SD and the CCSD methods recover all the diagrams for valence 
energies through the second order of MBPT, but start missing diagrams associated with valence triple excitations 
in the third order |8|]. Similarly, for contributions to matrix element of a one-body (e.g., electric dipole ) operator, 
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the SD method subsumes all the diagrams through the third order but misses approximately half of the diagrams in 
the fourth order of MBPT. The omitted fourth-order diagrams are entirely due to triple and disconnected quadruple 
excitations [24]] ■ Our group has carried out calculations of these 1,648 complementary diagrams for Na 25\ and 
Cs |17| . Close examination of our computed complementary diagrams reveals a high ( a factor of a hundred ) degree 
of cancelation between different contributions. Such cancelations could lead to a poor convergence of the MBPT series. 
Poor convergence calls for an all-order summation scheme and this is what we address here. The resulting formalism 
will recover the dominant fourth-order contributions to matrix elements and all third-order MBPT contributions to 
the valence energies in a nonperturbative fashion. 

The next systematic step in improving the SD method would be an additional inclusion of triple excitations 

Tc = 12 51 Pmnrabc <4n a t, a t a c a ba a , (14) 
mnrabc 

T-i) g ^ ^ Pmnrvab ^\ n Oj\ l 0^ r (li > CL a CL v (1^) 

mnrab 

into the cluster operator K (see Fig.0. However, considering the present state of available computational power, the 
full incorporation of triples (specifically, core triples) seems to be yet not practical for heavy atoms. 




FIG. 1: Diagrammatic representation of valence triple excitations. Double-headed arrow represents valence state. 



To motivate more accurate, yet practical extension of the SD method, we consider numerical results for the reduced 
electric-dipole matrix elements of 3s!/ 2 — 3pi/2 transition in Na 25] . From Table 1 of that paper, we observe that the 
contributions from valence triples T v (total —4.4 x 10~ 3 ) and non-linear doubles (disconnected quadruples) D n i (total 
1.3 x 10~ 3 ) are much larger than those from core triples T c (total 8 x 10~ 5 ). Similar conclusion can be drawn from 
our calculations for heavier Cs atom |17| . Because of this observation we will discard core triples and incorporate the 
valence triples into the SD formalism. We will refer to this method as SDvT approximation. Contributions of core 
triples to matrix elements are treated in this work perturbatively. 

In addition to triples, we will include effects from disconnected excitations. The relevant diagrams contribute at 
the same level as the valence triples and the full treatment of disconnected excitations will recover a part of the 
otherwise missing sequence of random- phase- approximation diagrams (see also discussion in Ref. ^3)- The resulting 
approximation will be referred to as CCSDvT method. 



III. FORMALISM 



Below we write down the CC equations for cluster amplitudes p in the CCSDvT approximation. The equations in 
the SD approximation are presented in Ref. |8j . We retain convention for the single and doubles from that paper and 
focus on additional terms due to valence triples and disconnected excitations. Some of the equations involving triple 
excitations were given in Ref. [l2l Il3j|: we use a different convention for the triples amplitudes. 



A. Valence triples 



In the following, we employ fully antisymmetrized valence triples amplitude p mn rvab- The object p m nrvab is anti- 
symmetric with respect to any permutation of the indices mnr or ab, e.g., 



'mnrvab Fnnirvab 



Pmnrvba Pmrnvba .... (1^0 
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It is straightforward to demonstrate that the contribution to the wave operator (and therefore all the resulting 
equations) can be expressed in terms of this antisymmctrized object. Explicitly, 

Tv — — ^ ^ Pmnrvab J \ n (l] fl (x! r (Xi ) (X a Cl v . (1 7) 

mnrab 

Computationally the use of pmnrvab substantially reduces storage requirements, as it is sufficient to store ordered 
amplitudes with m > n > r and a > b only. In the equations below, we will also use antisymmetrized combinations 
for doubles p mna b = Pmnab - Pmnba = Pmnab - Pnmab, Pmnva = Pmnva - Pnmva, and for the Coulomb matrix elements 

Qijkl — Qijkl Qijlk- 

From the general Eq.@ we obtain symbolically 

(£as ~f~ &b ~f~~ &v £m &n &r "~f~ $Ev) pmnrvab — 

T v [D c ] + T v [D v ] + T v [T v ] + T v [T c ] + nonlinear . (18) 

Here contribution T V [D C ] denotes effect of core doubles on valence triples, the remaining terms defined in a similar 
fashion. In this work we include only contributions T V [D C ] and T V [D V ] (see Fig. |3J) and omit the effect of valence 
and core triples on valence triples ( T^pl,] and T V [T C ]) and nonlinear CC contributions. Compared to the T V [D V ] and 
T V [D C ] contributions, these are higher order (and computationally expensive) effects. Explicitly, 

T V [D C ] 

^ (fJrncvaPnrcb QmcvbPnrca ~t~ QncvaPrmcb QncvbPrmca H - QrcvaPmncb QrcvbPmnca) 

c 

H~ ^ (fJnrsvPmsab ~r~ Qrmsv Pnsab fjmnsv Prsab ) , (19) 

S 

Tv^Dy] — ^ (jJrncabPrirvc ~i~ QncabPrmvc QrcabPmnvc) 
c 

^ {iJnrsbPrnsva QnrsaPmsvb ~r~ QrmsbPnsva QrmsaPnsvb QmnsbPrsva QmnsaPrsvb) ■ 

(20) 

s 

Notice that the matching of diagrams in Eq.© is generally not unique; we require that the r.h.s. of the above 
equation is fully antisymmetrized as the amplitude Pmnrvab on the l.h.s.; such procedure is unique and corresponds to 
a projecton of the CC equations onto the many-body state a^a^alabda^c) ■ Also from these equations we immediately 
observe that the triples enter the many-body wave function in the effective second order of MBPT, as the doubles 
enter in the first order in G, Eq. l|13|) . 



T V [D C ] T V [D V ] I 

FIG. 2: Representative contributions to the r.h.s of the valence triples equation. Horizontal dashed line denotes Coulomb 
interaction and solid lines — cluster amplitudes. 



B. Modifications to SD equations and valence energies 



Here we present CC equations for correlation energy SE V , valence singles p mv , and for valence double p m nva cluster 
amplitudes. In formulas below we write SD to denote contributions in the singles-doubles approximations tabulated 
in Refs. @>^2- As to the core amplitudes, they will be determined in the SD approximation (i.e. we do not include 
non-linear CC terms and core triples). 

Topological structure of the valence singles equation is 

\£v — £m + SE V ) p rnv = SD+ 

S v [S c <g> S v ] + S v [S c ® S B ] + S v [S c ® D v ] + S v [S v ® D c ] + S v [T v ] , (21) 
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where notation (K) n \{K) P ® (K) m ] stands for a contribution from a disconnected (p + m)-fold excitation (resulting 
from a product of clusters (K) p and (K) m ) to the cluster (K) n . We do not include cubic non-linear term S V [S C £g> 
S c ® S v ]. Explicitly, 



Sy [Sc ® £*v\ ^ ^j jaranr PnaPrv 7 

anr 

S V [S C ®S C ] = 

abnvPinaPnb ? 

abn 

Sy [Sc ® — ^ ffafenr (pmbPnrva PnbPmrva) 

abnr 

Sv [Sv ® -^c] ' ^ ^ Q abnr Pnv Pmrab j 



■ibirr 



Sv[T v ] = \^9a 

Representative diagrams are shown in Fig.[3| 



bnr Prnnrvab • 



abnr 



(22) 
(23) 
(24) 
(25) 
(26) 



S v [S c S v ] 



S V [S C ® S c 




FIG. 3: Sample contributions of triples and disconnected excitations to the valence singles equation. 



Valence doubles equation for p mn va can be symbolically represented as (see Fig. 

(e v + £ a — £ m — £ n + 5E V ) p m nva = SD+ 

D v [S c <gs S„] + D v [S c <g> Sc] + 

A, [5 C ® A,] + A [5„ ® A] + D v [S c ® A] + (27) 
A [A ® A] + A [5 C ® T u ] + A [S v ® T c ] + A [T„] . 

Contribution A [A ® £) c ] is topologically impossible and we omit cubic and higher-degree nonlinear terms like 
D V [S C ®S C ® S v ], D V [S C ® S c ® A], and D V [S V <g> S c <8 S c ® S^]. 



Explicitly, 



D V [D C ® D v ] = ^2,gbcrs ^yPrsvaP 



'mnbc "F ^Prasvapnrbc 



1- 

i~~p8nvaPmrbc PrsvbPnrnac T" PrsabPrnnvc 

~~ / j QbcTsPnhsvbPriTaci 
bcrs 



&v [$v & ^c] — ^ ^ Qbmrs Prv Pnsab H~~ ^ ^ Qbcar Prv Pnmbc j 
6rs 6cr 

[*^c & — 2 ^ ^ Qbnrs Prb Pmsva ^ ^ ^ fjbmrs Prb Pnsva 

brs brs 

H~ 2 / ^ Qbmrs Pnb Prsva ^ / J 9brirs Prnb Prsva 
brs brs 

^ Qbnrs Pra Pmsvb / Qbcar Pre Pmnvb ^ Qbcar Pnc Prrnvb ) 
6rs bcr bcr 

&v [$c ® ^c] — ^ ^ Qbcvr Pre Pnmab ^ ^ Qbcvr Prnc Prnab ~t~ ^ ^ Qbcvr Pra Prnnbc j 
6cr 6cr bcr 

[5*c ® S v \ — ^ Qbnar Pmb Prv "t - ^ ^ Qmnrs Prv Psa ? 
br rs 

D V [S C ® S c ] =J2 

Qbmvr Pnb Pra \ / Qbcav Prnc Pnb • 
br be 

The effect of valence triples on valence doubles reads 

~~ 2 ^ ^ {ffbcarPmnrvbc QbcvrPnmrabc) 
rbc 

+ (9bnrs Pmsrvab \ QbrnrsPsnrvab) ■ 

rsb 




D V [D C ®D V ] a 

FIG. 4: Effects of disconnected and valence triple excitations on valence doubles. 

Finally, the valence correlation energy may be represent as 

SE V = 5E SB + 5E CC + 5E vT , (28) 

with 

SE CC = 

Qavnr PnaPrv \ / QabnvPvaPnb 

(29) 

anr abn 

"F ^ ^ Qabnr [PvbPnrva PnbPvrva Pnv Pvrab] , 
abnr 

SE V T — — gabmnPvmnvab ■ (30) 



abn in 



Topological structure of contributions to energy is similar to the terms on the r.h.s of the valence singles equation l|21|l . 
Here correction 8Eqc comes from non-linear CC contributions and 6E v t is due to valence triples. 
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C. Normalization 

The CC wave function is derived using the intermediate normalization, \^v) — 1 and in calculating the atomic 
properties based on the CC wave function, one needs to renormalize it. In calculations of matrix elements one requires 
the valence part of the normalization, N v — (4 r „|W„)vai, connected- We obtain 

— SD -|- ^ ^ Pmnab Pvmnvab ~\~ ^ ^ {.Pmnrvab) ■ (31) 

mnab mnrab 

The last term in the equation above is quadratic in valence triples (i.e., it is of the fourth effective order) and we will 
neglect it in the following. 



D. Matrix elements of one-body operator 

Finally, we consider matrix elements of a one-body operator Z = Ylij z ij a \ a i between two CC states \*&v) and 
1^^,). Taking into account renormalization, this matrix element can be defined as 

_ (9 W \Z\9 V ) 

As it was shown in Ref. |8( all disconnected diagrams in the numerator and denominator of this expression cancel, 
leading to 

M wv = V wv,conn ttt • (33) 

{[l + ^Unn] [l + OWcoJ) ' 

We discarded valence-independent contribution, as it vanishes for non-scalar operators. To unclutter the notation 
below we simply write 

Zwv — (^u)comi) 
N v = (K al )conn. (34) 

Blundell et al. Q tabulated 21 contributions to the matrix elements in the SD approximation. These SD corrections 
are mainly due to (i) random-phase-approximation (RPA) diagram proportional to a product of Z and D v and (ii) 
the Brueckner-type (core-polarization) diagram proportional to the product of Z and S v . In Ref. |17| we additionally 
included modifications to M. wv caused by non-linear terms in the CC wave function. We have devised a re-summation 
scheme that is equivalent to "dressing" of lines and vertices of the SD diagrams (see also Ref. p6j). 

Including valence triples leads to additional direct contributions, Z wv = SD + Zwv ■ We obtain 

7 

Zwifi = Z^v ,k ^ ■> (35) 
fc=i 

Zwv ' ^ — ^ PmaPwmnvab z bn ~t~ h.C.S., (36) 

abmn 

Zwv' ^ = ~2 ^ y PmnbaPwnmvcb z ca ~\~ h.C.S., (37) 

abemn 

Z wv v = — ^ ] P m nabPrmnvab z wr ^~ h.C.S., (38) 
abmnr 

Z W y' ' = — ^ ] P m nabPwrmvba z nr + h.C.S., (39) 
abmnr 

Z W y' ' = — — ^ t PranwbPrmnvab z ar + h.C.S., (40) 
abmnr 

Zyjv ' = — r y ] PmnrwcbP~mnrvab z ac, (41) 
abemnr 

Z W v" — 7 ^ y PmnrwabP 'snrvab z ms ■ (42) 



4 

abmnr s 
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In these expressions, abbreviation h.c.s. stands for a Hermitian conjugation of the preceding term with a simultaneous 
swap of the valence indices w <-> v. As discussed in Ref. |24|. valence triples start contributing in the fourth order of 

(Tv k) 

MBPT for matrix elements; these contributions correspond to terms Z wv ' , k = 2 — 5. We presently discard terms 
=#=6 and #7 that are quadratic in triple excitations. 



E. Symmetries and reduced triples 



Relativistic one-particle orbitals i are characterized by the principle quantum number rn , the total angular momen- 
tum ji, it's projection rrii and the orbital angular momentum U. The summations over magnetic quantum numbers 
are carried out analytically substantially reducing the number of coefficients. A dependence of valence triples on 
magnetic quantum numbers may be parameterized as (we use angular diagrams, see, e.g., Ref. |23jp 



J mnrvab 



E 

LL'h 



]nm n 



■j v m v 



j,.m r 



\ 1 

L ) 


* V ' 


>— 

+ 


5» 



+ 



FhL'h (mnrvab) 



(43) 



Km 



where h is a half-integer coupling angular momentum and L and L' are integer coupling momenta. The "reduced 
triples" Fll' h (mnrvab) do not depend on magnetic quantum numbers. 

Selection rules for various angular momenta characterizing reduced triples follow from properties of the 3j-symbols 
in the angular diagram (|43[) . In addition, the atomic Hamiltonian is invariant under parity transformation, leading 
to an additional parity selection rule l m + /„ + l r + l v + l a + 4 = even integer for a triple amplitude p m nrvab ■ 

Owing to the antisymmetric properties of the triples, Eq. l|16ll . it is sufficient to store reduced triples with (n m x m ) > 
(n n K n ) > (n r K r ) and (n a x a ) > (n^xt,), where a = (I — j)(2j + 1). The reduced triples with other combinations of 
arguments can be related to the ordered set via symmetry properties. For example, 

( jb h L) 

F LL >h (mnrvba) = (2h + 1) (22/ + 1) ^ I j r L' j a \ (-l) h+h +K+L F LKh , (mnr vab) . (44) 

h'K \ K j n hi J 

There are 11 such index-swapping relations for reduced valence triples. 



IV. NUMERICAL RESULTS AND DISCUSSION 



To reiterate discussion so far, we derived algebraic expressions in the CCSDvT formalism, which includes valence 
triples and a subset of disconnected excitations. We also carried out angular reduction of these expressions and 
developed a numerical code. In this section we present our ab initio results for properties of 3s, 3p!/ 2 , and 3p 3 / 2 
states of atomic sodium. Results for removal energies are presented in Section IIV Al and for dipole matrix elements 
and HFS constants A in Section TlV Bl 

Before presenting the results, let us briefly describe our numerical code. It is an extension of the relativistic SD 
code [T^| which employs B-spline basis set. This basis numerically approximates complete set of single-particle atomic 
states. Here we use 35 out of 40 positive-energy (s.i > — m e c 2 ) basis functions. Basis functions with Z max < 6 are used 
for singles and doubles. For triples we employ a more limited set of basis functions with l max (T v ) < 4. Excitations 
from all core sub-shells are included in the calculations. Numerically we found that this choice is a reasonable 
trade-off between storage and overall numerical accuracy (after all, triples affect computed properties at ~ 1% level.) 
The results presented in this Section will include basis set extrapolation correction, which is obtained by computing 
SD properties with increasingly larger basis sets and interpolating them to I = oo. The CC equations were solved 
iteratively. We notice that the reported calculations can be carried out in the memory of a modern high-end personal 
workstation: storing reduced valence triples in a single precision required about 900 Mb for Sx/a states and 1.5 Gb 
for 2? 3 /2 states (the latter involve more angular channels). 
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A. Energies 

Computed removal energies of 3s, 3pi/2, and 3^3/2 states of atomic sodium are presented in Table [I] The dominant 
contribution to the energies comes from the DHF values. The remaining (correlation) contribution is given by Eq. 
We computed this correlation correction in several approximations: SD, SDvT, CCSD, and, finally, CCSDvT. 



TABLE I: Contributions to removal energies of 3s, 3pi/2, and 3p3/2 states for Na in cm in various approximations. A 
comparison with previous CC-type calculations and experimental values is presented in the lower panel. 





3s 




3pi/2 


3P3/2 


-Edhf 


39951.6 




24030.4 


24014.1 


SE so 

TTltot 


1488.8 
41440.3 


SD 


463.9 
24494.3 


460.6 
24474.7 


5E^ iT 
5E yT 

TTttot 

-^SDvT 


79.7 
25.4 
41545.5 


SDvT 


28.9 
4.8 
24528.0 


28.4 
4.7 
24507.8 


SEcc 

TTttot 

-C'CCSD 


-57.0 
-17.5 
41365.9 


CCSD 


-20.0 
-7.4 
24466.9 


-18.4 
-7.4 
24448.9 


SE vT 
SEcc 

TTttot 

-C'CCSDvT 


16.8 
23.7 
-18.4 
41462.5 


CCSDvT 


6.8 
4.5 
-8.0 
24497.6 


7.9 
4.4 
-8.0 
24479.1 


SD(pvT) [13] 
CCSD [iOj 

rp a 
J-^cxperim 


41447.3 

41352 

41449.6 


Other works 


24493.9 

24465 

24493.4 


24476.7 
24476.2 



These values are from spectroscopic data compiled bv NIST |27l| 



First we list correlation energies SEg^) obtained in the SD approximation. The results contain basis set extrapolation 
corrections from Ref. |2^. The extrapolation corrections increase the removal energies by 5.1 cm" 1 for the 3s state, 1.9 
cm" 1 for the 2>P\/2 state, and 0.8 cm -1 for the 3p 3 / 2 - Total removal energy is = -Edhf +<5-Esd- At the next step 
(SDvT) we include valence triple excitations, i.e., in the CC equations in addition to the SD terms we incorporate terms 
with amplitudes Pmnrvab- It is instructive to distinguish direct and indirect SE 1 ^" effects of these excitations. The 
direct effect of triples is 5E v t, Ea. (|30|l . while indirect effect is a modification of SE^b due to effect of triples through 
coupling to singles and doubles. In this case, the indirect contribution is defined as SE^" = SEsu [SDvT] — SEst> [SD] . 
We list the two types of contributions in the Table and it is clear that for all the states both contributions add 
constructively, and for all the considered approximations the indirect contribution dominates over the direct one. The 
total removal energy in the SDvT approximation is #sdvT = e dhf + SE SB + SE^ ir + SE vT . The totals for other 
approximations are defined in a similar way. 

As we move to the CCSD approximation in Table [I] we notice that here the corrective terms 8E<^ lr and 5Eqq 
decrease the removal energies, while for the SDvT case the corrections increased E tot . In both cases the resulting total 
energies E tot were moved away from the experimental values. Since the effects of disconnected and triple excitations 
are comparable and opposite in sign, simultaneous treatment of the two effects is required. The results of such 
treatment are listed under CCSDvT heading in the Table. Compared to the CCSD and SDvT approximations, the 
CCSDvT results move into a closer, 0.01-0.03%, agreement with the experimental values. 

Comparison with the previous CC-type calculations of Na removal energies is presented in the lower panel of the 
Table [I] SD(pvT) approximation denotes results obtained with a scheme originally proposed in Ref. [jj. In this 
scheme: (i) starting from the SDvT approximation, one keeps vT contributions in the equation for valence singles 
and valence energies (i.e., D V [T V ] effect is neglected) ; (ii) triples are approximated by 
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(iii) to avoid expensive storing of valence triples, in the p mv equation the triples denominators 
(e a + £5 + e v — e m — e n — e r + SE V ) are replaced by an approximate combination (e a + — e„ — e r ). In this ap- 
proximation S„[T„] effect is effectively overemphasized ( for the ground state e v < e m ). In the expression for the 
energy, 5E v t, Eq.J3U|), triples enter as p V mnvab and the above replacement of denominators is more algebraically 
justified. Nevertheless, we found a substantial (a factor of three) disagreement between SE v t corrections obtained in 
our (more complete) SDvT and SD(pvT) approximations. 

To understand the origin of this large disagreement, we have compared individual contributions to 5E v t coming 
from the r.h.s. of the triples equations with the corresponding contributions in the SD(pvT) approximation. We found 
that the individual terms agree at a reasonable 10% level. The discrepancy in the total value arises because there 
are certain very large individual terms canceling each other. These terms are several hundred times larger then the 
final combined result. In other words there is a subtle cancelation taking place and our more sophisticated all-order 
treatment profoundly affects this delicate cancelation. 

In addition, in Ref. [l^. the explicit contributions of triples to the energies, 6E v t, were computed using direct 

third-order MBPT approach. Such terms are denoted in Ref. 12] as -E^cxtrai to emphasize that these are diagrams 

(3) 

missed in the SD approximation in the third order. A comparison of our computed SE V ^ with E^ <3 Xtra is presented in 

Table [H] We again observe a large discrepancy, due to substantial cancelations among contributions to £^ xtra and 
resulting enhanced sensitivity to a correct all-order treatment. 

The CCSD results obtained by Eliav et al. [13 agrees with our CCSD energies for the Spi/2 state. However, for 
the 3si/2 the two calculations disagree by 14 cm" 1 . This discrepancy is likely due to our omission of all non-linear 
terms in the core CCSD equations. 



TABLE II: Comparison of complementary third-order MBPT contributions E v oxtra to removal energies with the corresponding 
all-order correction SE v t- The corrections are given in cm -1 . 





3s 


3fl/2 


3^3/2 


8E vT 


-25.4 


-4.8 


-4.7 


Ei% ia , Ref. [12] 


-9.2 


-1.5 


-1.6 



Comparing the final CCSDvT results for the removal energies with the experimental values (last row of Table |nj) 
we find an agreement at the level of O.OTO.03%. We do not include Breit-, reduced-mass and mass-polarization 
corrections to the energies, as they contribute at a much smaller level M A perfect theory-experiment agreement for 
the previous SD (pv T) calculations of energies is fortuitous because contributions of the disconnected excitations 
omitted in Ref. [2!j would move the theoretical energies by about 70 cm -1 for the 3s!/ 2 state (see Table [IJ|. 



B. Hyperfine constants and electric-dipole amplitudes 

With the computed wave functions of the 3s, ?>Pi/2 and 3p 3 / 2 states we proceed to determining magnetic-dipole 
hyperfine-structure constants A and electric-dipole transition amplitudes. The formalism was outlined in Sec. IIII Dl 
and here we discuss our ab initio results and compare them with the experimental values. 

Numerical results are presented in Table 11111 This Table is organized as follows. First we list the DHF and SD 
values. The results for the HFS constants include finite- nuclear size effects (see Appendix). In the part denoted 
"All-order corrections beyond SD" we tabulate differences between the values obtained at a certain approximation 
(CCSD, SDvT, CCSDvT) and the corresponding SD value (symbolically, e.g., A(CCSD) = CCSD - SD). The most 
sophisticated approximation is CCSDvT (it includes both implicit and explicit, Eq. 1351) . contributions of valence 
triples and implicit contribution of disconnected excitations); we will base our final ab initio result on the CCSDvT 
values. A cursory look at this part of the Table reveals that the contributions of disconnected excitations tend to 
compensate contributions of valence triples for all the computed properties. This situation is similar to the one 
observed by us while presenting results for removal energies in Section llV Al 

While discussing the CCSDvT results, it is instructive to compare the explicit valence triple corrections to matrix 
elements, Ea. ()35|l . with a corresponding contribution from the direct fourth-order calculations 25]. In particular, for 

the (3s||D||3pi/ 2 ) amplitude, the Z^v' CCSDvT contribution of -0.00075 is in a close agreement with the fourth-order 
Zix2(T v ) contribution of -0.00073. The close agrement is due to the fact that there are no strongly canceling terms 
in the Z\ X 2 {T v ) class of the fourth-order diagrams. This should be contrasted with our similar comparison of energy 
corrections (see Table ITTf) . where large, a factor of 100, cancelations lead to a poor accuracy of the direct third-order 
computation. 
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TABLE III: Hyperfine structure constants A (in MHz) and matrix elements of electric dipole moment (in a.u. ) for 23 Na. 
Results of calculations and comparison with experimental values. See text for the explanation of entries. 





A(3s) 


A(3p 1/2 ) 


^(3p 3 / 2 ) 




\ t) P3/2 1 1 — 1 


DHF 


623.8 


63.39 


12.59 




5.2188 


SD 


889.0 


95.05 


18.85 


3.5308 


4.9932 






All-order corrections beyond SD 






AfOOSD"! 


—7.7 


-1.76 


-0.34 


0.0072 


0098 

W . WW CO 


A(SDvT) 


8.6 


2.06 


0.36 


-0.0115 


-0.0166 


A(CCSDvT) 


0.4 


0.07 


-0.02 


-0.0035 


-0.0053 






Complementary corrections 






Line dressing 


-2.4 


-0.43 


-0.09 


0.0004 


0.0005 


Vertex dressing 


1.5 


0.17 


0.04 


-0.0001 


-0.0002 


MBPT-IV (core triples,...) 


-2.8 


-0.41 


-0.06 


0.0001 


0.0001 


Breit + QED [5, 30] 


0.2 






0.0001 


0.0002 


Final CCSDvT + corrections 


885.9 


94.45 


18.72 


3.5278 


4.9885 


Experiment 


885.81 a 


94.44(13) 6 


18.534(15) c 


3.5267(17) d 


4.9875(24) d 



94.42(19) e 18.572(24/ 3.5246(23) 9 4.9839(34) 9 

18.64(6)'' 
18.69(9)* 

Agreement with experiment 0.01% < 0.1% 1% < 0.05% < 0.05% 

"Ref. 31]; 6 Ref. 32]; c Ref. 18]; d Ref. 33]; e Ref. 34]; 'Ref. 19]; 9 Rcf. 35]; A Rcf. 20]; ! Rcf. 21]. 

Corrections beyond the CCSDvT approximation are listed in the Table IIIII under the heading "Complementary 
corrections". The dressing corrections arise due to a direct contribution of disconnected excitations to the matrix 
elements. The details of our all-order dressing scheme can be found in Rcf. Following that work we distinguish 
between vertex- and line- dressing corrections. Futher, the "MBPT-IV" entries in the Table include all IVth diagrams 
missed by the CCSDvT method and dressing. For example, our CCSDvT approximation discards core triples and 
disconnected core excitations and these contributions arise starting from the fourth order of MBPT for matrix elements. 
In notation of Ref. [24| the complementary fourth-order terms are Zq X 3(D v {T c ]), Zq X 3(S c [T c ]), and Zi X 2(T c )- In 
addition, the dressing method of Ref.[l7| misses so-called stretched and ladder Z\ x 2 {D n i) diagrams. These diagrams 
are also incorporated into the value of "MBPT-IV" contribution in the Table ITTT1 We used the fourth-order code of 
Ref. (2^] to evaluate the complementary MBPT-IV contributions. 

Finally, we tabulate Breit and QED corrections available from the literature (see Appendix for discussion). By 
combining them with the CCSDvT values and the complementary corrections we arrive at the final ab initio values in 
the bottom part of Table ITTT1 Here we also present a comparison with the experimental data. In particular, the last 
row tabulates percentage deviations from the experimental values. If the ab initio value lays inside the experimental 
error bar, we tabulate experimental uncertainty instead. The theory-experiment agreement is better than 0.1% except 
for the HFS constant of the 3^3/2 state, where our value disagrees with most accurate experimental results at 1% 
level. For this constant our result is, however, in a reasonable agreement with the less accurate (0.3% uncertainty) 
result of Ref. 

V. CONCLUSION 

To reiterate here we presented a practical high-accuracy ab initio relativistic technique for calculating properties of 
univalent atomic systems. The distinct formal improvements over the previous singles-doubles approach |8l M IT2I Il3| 
are: 

1. non-perturbative treatment of valence triple excitations; 

2. incorporation of disconnected excitations (non-linear terms) in the coupled-cluster approach; 

3. inclusion of complementary MBPT diagrams so that the calculations of matrix elements are complete through 
the fourth-order of MBPT; these diagrams include contributions of core triples. 

4. all-order "dressing" of lines and vertices in calculations of matrix elements. 

Including all the enumerated effects is important in reaching the present uniform "better than 0.1%" theoretical 
accuracy for Na atom. In particular, a simultaneous treatment of triple and disconnected quadruple excitations is 
required, as these two relatively large effects tend to partially cancel each other. 



13 



In the framework of the developed formalism, we computed removal energies, magnetic-dipole HFS constants A and 
electric-dipole amplitudes for the principal 2>Si/ 2 — Zpj transitions. The presented approach demonstrates a uniform 
sub-0.1%-accurate agreement with experimental data. In particular, we find that the removal energies are reproduced 
within 0.01-0.03% and the HFS constants of the 3s and ip\/2 states with a better than 0.1% accuracy. The calculated 
dipole amplitudes are in a perfect agreement with the 0.05%-accurate experimental data. In the case of the 3p 3 / 2 state 
HFS constant our ab initio result deviates from ~ 0.1%-accurate experimental values [UBil by 1%, while agreeing 
with the less accurate measurements |20t \'2\\ . We anticipate that the relativistic many-body techniquepresented here 
can serve as a basis of highly-accurate evaluation of parity- violating effects in Cs atom and Ba + ion [4( . 
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APPENDIX A: SMALLER (NON-CORRELATION) CORRECTIONS TO THE HYPERFINE 

STRUCTURE CONSTANTS 

Calculations of magnetic hyperfine constants A presented in Tablc lTTTl were carried out with the nuclear gyromagnetic 
ratio gi — 1.4784. In calculations we model the nucleus as a uniformly magnetized sphere of radius 3.83 fm. For the 
3si/ 2 state, the corresponding nuclear size (Breit-Weisskopf) effect reduces point-nucleus results by 0.5 MHz. In an 
extreme case, when magnetization is assumed to be completely localized on the nuclear surface, the Ahf s (3si/ 2 ) is 
further reduced by 0.15 MHz; this difference between the uniform and surface magnetization is below our theoretical 
accuracy. 

Breit and QED contributions to the HFS constant of the 3si/ 2 state were calculated recently by Sapirstein and 
Cheng |{|. In their notation, the value marked "Breit/QED" includes effects of the Breit interaction, retardation in the 
transverse photon exchange and negative-energy states, while "QED" correction encapsulates vacuum polarization 
and self-energy corrections. (The Breit correction of 0.35 MHz, evaluated using analytical expression |36( is in a 
reasonable agreement with the value of 0.2 MHz from Q). As to the QED corrections, the leading Schwinger term 
(anomalous magnetic moment) SA/A = a/n sets a scale for radiative corrections at 0.1% and this is comparable with 
the accuracy of our calculations. Nevertheless, explicit model-potential calculation |5j of vacuum polarization and 
self-energy corrections displays a large degree of cancelation between different contributions, leading to the total QED 
correction 70 times smaller than the Schwinger term. 

Following discussion of Ref. [37| for Li, we also analyzed the following smaller corrections to the HFS constant: (i) 
Mass scaling. This effect contributes at the relative level of 1/(1 + m e /M nuc ) 3 ps 7 x 10~ 3 %; here M nuc is the nuclear 
mass, (ii) Mass polarization. It occurs due to an additional introduction of the term — /i/M nuc J2i>j ^« ' ^ j the 
atomic Hamiltonian, [i being the reduced mass of the electron. We expect that this term would contribute at the 
relative level of l/M nuc (aZ) 2 10 -5 %. (iii) Second order in magnetic-dipole HFS interaction. It contributes at the 
10 -5 % level. All the enumerated corrections are below the level of theoretical accuracy of the calculation presented 
here. 
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